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ABSTRACT 

We  construct  a  free  energy  density  functional  approximation  for  the  primitive  model 
of  the  electrical  double  layer.  The  hard  sphere  term  of  the  free  energy  functional  is  based 
on  a  nonlocal  generic  model  functional  proposed  by  Percus.  This  latter  model  functional, 
which  is  a  generalization  of  the  exact  solution  for  the  non-uniform  hard  rod  model,  requires 
as  imput  the  free  energy  of  a  homogeneous  hard-sphere  mixture.  We  choose  the  extension 
of  the  Carnahan- Starling  equation  of  state  to  mixtures.  The  electrostatic  part  of  the 
non-uniform  fluid  ion-ion  correlations  present  in  the  interface,  is  approximated  by  that  of 
sin  homogeneous  bulk  electrolyte.  Using  the  mean  spherical  approximation  for  a  neutral 
electrolyte,  we  apply  the  theory  to  symmetrical  1:1  and  2:2  salts  in  the  restricted  primitive 
model.  We  present  comparisons  of  density  profiles  and  diffuse  layer  potentials  with  Gouy- 
Chapman  theory  and  Monte  Carlo  data.  When  available,  we  also  compare  our  results 
with  data  from  other  recent  theories  of  the  double  layer.  For  highly  charged  surfaces,  the 
profiles  show  the  layering  of  counterions  and  charge  inversion  effects,  in  agreement  with 
Monte  Carlo  data. 
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I.  INTRODUCTION 


Understanding  the  behavior  of  charged  particles  near  charged  surfaces  is  an  impor¬ 
tant  problem  in  physical  chemistry.  Separation  of  charge  in  response  to  the  field  of  the 
charged  surface  is  refered  to  as  the  electrical  double  layer.  Double  layers  are  present  in 
electrochemistry  in  the  form  of  the  electrode/electrolyte  interface,  and  they  often  play  a 
major  role  in  the  stability  of  soap  films,  colloidal  dispersions,  and  biological  membranes. 
As  a  result  of  the  occurrence  of  double  layers  in  numerous  situations,  there  has  been  a  con¬ 
siderable  effort  to  describe  them  theoretically.  The  early  theory  that  met  with  significant 
success  was  that  of  Gouy  [1]  and  Chapman  [2]  and  was  based  on  the  Poisson- Boltzmann 
equation.  More  recently,  theory  has  been  built  on  more  rigorous  methods  of  the  statistical 
mechanics  of  the  liquid  state  [3].  Because  the  physical  systems  in  which  double  layers 
occur  are  generally  quite  complicated  (for  a  recent  review  see  Ref.  [4])  theoretical  efforts 
have  been  directed  towards  determination  of  the  properties  of  greatly  simplified  models. 
The  Gouy-Chapman  theory  [1-2]  was  developed  for  a  model  of  point  charges  next  to  a  uni¬ 
formly  charged  planar  surface,  for  example.  A  later  modification  of  this  theory  by  Stem 
[5]  ,  known  as  the  modified  Gouy-Chapman  theory  (MGC),  is  based  on  the  same  model. 
This  description  is  quite  accurate  provided  the  real  ionic  radius  is  not  too  large  compared 
to  ionic  spacing  and  the  charge  on  the  surface  and  on  the  particles  are  relatively  small 
(low  density-weakly  coupled  systems).  At  higher  densities  or  for  highly  coupled  systems, 
the  core  interaction  becomes  important.  Thus,  to  take  into  account  of  the  finite  size  of 
the  charged  particles,  many  authors  focused  their  attention  on  a  model  electrical  double 
layer  composed  of  charged  hard  spheres  at  a  hard,  planar,  polarizable,  uniformly  charged 
surface.  The  ions  are  assumed  to  be  immersed  in  a  continuum  with  a  dielectric  constant 
which  may  be  different  from  that  of  the  charged  wall.  This  model  is  known  as  the  primitive 
model  (PM)  of  the  double  layer. 


There  is  a  considerable  body  of  recent  work  on  the  PM  double  layer.  For  testing 
theory,  the  Monte  Carlo  (MC)  simulations  by  Valleau  and  co-workers  [6],  are  especially 
significant.  There  is  the  work  on  the  modified  Poisson- Boltzmann  (MPB)  approximation 
[7],  which  is  based  on  the  Kirkwood  [3]  hierarchy.  In  that  approach  the  effects  of  the 
wall  on  the  ion-ion  correlations  are  handled  in  a  natural  way.  On  the  other  hand,  the 
work  based  on  the  singlet  Ornstein-Zernike  (OZ)  equation  [8]  with  the  mean  spherical 
approximation  (MSA)  or  the  hypemetted  chain  approximation  (HNC)  stresses  on  a  more 
careful  treatment  of  the  effects  due  to  the  finite  size  of  the  ions  while  the  direct  ion-ion 
correlation  functions  near  the  surface  axe  approximated  by  the  functions  calculated  in  the 
bulk  solution  [9].  There  is  also  work  based  on  the  Born- Green-Y von  (BGY)  equation 
[10-11].  By  using  phenomenological  expressions  for  the  inhomogeneous  pair  correlation 
functions  as  closures  for  the  BGY  equation,  this  latter  method  emphasizes  the  importance 
of  properly  handling  the  ion-ion  correlation  functions  near  the  wall.  The  BGY  equation 
exactly  satisfies  the  contact  theorem  [3].  This  is  especially  important  in  the  high  density- 
high  coupling  regime.  The  recent  work  by  Fortsmann  and  collaborators  [12,13],  in  which 
the  ion-ion  direct  correlation  functions  are  computed  using  the  MSA  at  a  local  non-neutral 
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concentration  (HNC/LMSA),  is  also  aimed  at  building  into  the  theory  good  pair  correlation 
functions.  Similar  in  spirit  to  the  work  of  Fortsmann  and  collaborators  is  that  of  Kjellander 
and  Marcelja  [14],  in  which  the  double  layer  interaction  between  two  uniformly  charged 
surfaces  inmersed  in  an  electrolyte  solution  is  calculated.  Perhaps  the  most  accurate  recent 
work  is  that  of  Plischke  and  Henderson  [15].  In  that  work,  the  inhomogeneous  OZ  equation 
with  the  HNC  and  MSA  closures  were  solved  together  with  the  Lovett-Mou-Buff-Wertheim 
equation  [16]  for  the  density  profiles  of  the  ions  (OZ/LMBW). 


Double  layers  are  good  examples  of  strongly  inhomogeneous  systems.  In  an  elec¬ 
trolyte,  especially  at  elevated  surface  charges,  the  density  variations  near  the  electrode  are 
extremely  large.  The  situation  of  an  ion  near  the  wall  is  totally  different  from  that  of  a 
similar  ion  in  the  neutral  bulk.  Ideally,  a  double  layer  theory  should  take  into  account  the 
correlations  arising  from  both  the  hard-core  repulsion  and  the  electrostatic  interactions. 
Because  these  correlations  are  strongly  dependent  on  the  distance  from  the  wall,  only  a  few 
theories  are  able  to  handle  properly  the  ion  packing  near  the  electrode.  According  to  MC 
results  for  1:1  electrolytes  [6], at  high  electrode  charges  the  counterions  start  the  formation 
of  a  second  layer  before  the  first  layer  is  densely  packed.  Of  the  theories  mentioned  above, 
only  the  HNC/LMSA,  the  BGY,  the  Kellander  and  Marcelja  and  the  OZ/LMBW  theories 
are  able  to  predict  the  formation  of  the  second  layer  of  counterions. 

Parallel  to  the  development  of  the  double  layer  theories,  the  last  decade  has  seen  a 
geat  deal  of  activity  in  the  study  of  non-uniform  fluids  using  free  energy  density  functional 
theories.  This  method,  which  originated  with  van  der  Waals  [17],  requires  the  construction 
of  an  expression  for  the  free  energy  of  the  inhomogeneous  system.  Even  though  the  rigor¬ 
ous  statistical  mechanics  formalism  of  density  functional  theory  was  established  more  than 
twenty  years  ago  [18,19],  the  reduction  of  the  exact  results  to  tractable  accurate  approxi¬ 
mations  has  been  the  goal  of  many,  more  recent  investigations  [20-26].  Treatments  based 
on  local  density  approximations  have  proven  useful  to  describe  weakly  structured  systems, 
like  fluid- fluid  interfaces  [20-22],  or  fluids  in  weak  external  fields,  but  are  not  applicable  to 
the  strong  inhomogeneities  characteristic  of  fluid-solid  interfaces.  In  order  to  handle  the 
strong  inhomogeneities  present  in  fluid-solid  interfaces,  a  nonlocal  approach  was  introduced 
by  Nordholm  and  coworkers  [23]  in  their  generalized  van  der  Waals  theory  (GVDW).  Since 
then,  systematic  improvements  in  the  method  in  which  finite  size  effects  are  considered 
have  been  published  [24,26].  Applications  of  nonlocal  theories  include  the  studies  of  the 
structure  of  confined  fluids  [27-30],  capillary  condensation  [31,32],  layering  transitions  [33], 
and  the  planar  electrical  double  layer  [34].  Very  recently,  the  nonlocal  smoothed-density 
approach  (SDA)  due  to  Tarazona  [24],  was  extended  to  binary  hard  sphere  mixtures  with 
arbitrary  size  ratio  [35]. 

In  this  work,  we  present  a  theory  for  the  electrical  double  layer  in  which  the  effects  due 
to  the  finite  size  of  the  particles  are  considered  within  the  framework  of  a  generic  nonlocal 
functional  proposed  by  Percus  [36]  and  generalized  to  multicomponent  fluids  by  Vanderlick 
et  al.  [29].  This  generic  functional  can  be  used  to  generate  the  functionals  of  several  known 
nonlocal  approaches  [30].  These  include  the  GVDW,  the  SDA,  and  a  functional  proposed 
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by  Robledo  and  Varea  [37],  and  by  Fischer  and  Heinbuch  [38],  as  a  generalization  of  Percus’ 
exact  solution  of  the  one  dimensional  hard-rod  system  [39].  We  use  the  latter,  termed  here 
generalized  hard-rod  model  (GHRM),  to  construct  a  model  functional  for  the  electrical 
interface.  Its  advantage  over  the  OZ/LMBW  theory  used  by  Plischke  and  Henderson  is 
that  its  working  equations  axe  much  cheaper  to  evaluate. 

The  article  is  as  follows.  The  PM  of  a  planar  double  layer  is  described  in  Sec.  II.  The 
general  free  energy  density  functional  formalism  for  the  electrical  double  layer  is  presented 
in  Sec.  III.  In  Sec.  IV  we  report  our  results  for  the  density  profiles  and  electrostatic 
potential  and  compare  these  with  the  MC  results  of  Valleau  and  co-workers  [6]  and  with 
results  of  some  of  the  theories  mentioned  above. 


II.  PRIMITIVE  MODEL 

In  the  primitive  model  of  the  electrical  double  layer,  the  electrolyte  is  assumed  to  be 
a  fluid  of  charged  hard  spheres  of  charge  qa  and  diameter  da  inmmersed  in  a  dielectric 
continuum  of  dielectric  constant  e.  Separating  the  Coulombic  and  short-range  repulsive 
contributions  to  the  pair  interaction,  we  have 


O  =  qaq&uc{\  r  -  r'  |)  +  uTa0(\  r  -  r'  |), 

(2.1) 

where 

uc(r)  =  1/er  , 

(2.2) 

and 

««*(»■)  =  r  <  (da  +  d0)/2 

=  0,  r  >  (da  +  dg)/2. 

(2.3) 

The  electrode  is  considered  to  be  an  infinite  flat  hard  wall  with  a  uniform  charge  density 
cr.  This  impenetrable  hard  wall  produces  a  repulsive  potential,  for  particles  of  species  a, 
of  the  form 

v;(x)  =  oc,  x<da/ 2, 

=  0,  x  >  da/ 2,  (“'4) 

where  x  is  the  ion’s  distance  to  the  plate.  On  the  other  hand,  the  uniform  surface  charge 
density  gives  rise  to  a  Coulombic  potential  of  the  following  form 

uc(x)  = -27r<7  |  x  | /e  +  C,  (2.5) 

where  C  is  a  constant  which  depends  on  the  choice  of  the  point  of  zero  potential. 

In  order  to  eliminate  image  charges,  the  dielectric  constant  is  also  taken  to  be  e  in  the 
region  x  <  da  /2.  The  total  external  potential  can  now  be  written  as 

Va(x)  =  qavc{x)  -(-  v;(x).  (2.6) 
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A  quantity  of  an  enormous  importance  in  the  electrical  double  layer  theory  is  the  mean 
electrostatic  potential.  The  mean  electrostatic  potential  t/>(r)  at  a  point  r  is  related  to  the 
density  distribution  functions  na(x)  in  the  following  way 

4 

0(r)  =  uc(x)  +  J  d3r'uc(\r -r'  \)^2qana(x').  (2.7) 

The  formal  solution  to  Poisson’s  equation  yields  the  following  expression  for  the  mean 
electrostatic  potential 


a 


(2.8) 


The  boundary  conditions  used  in  arriving  at  Eq.(2.8)  axe  0(oo)  =  0,  and 


di/>(x)  .  —4 iru 

— 1 —  r=o=  - 

ax  6 


In  the  derivation  of  Eq.  (2.8)  we  required 


roa 

/  dx'  Y^qc.na(x')  =  -tr, 
Jo 


(2.9) 


(2.10) 


which  is  the  constraint  of  overall  electroneutrality  of  the  system.  From  Eq.  (2.8)  we  can 
observe  why  the  mean  electrostatic  potential  evaluated  at  the  closest  approach  distance  is 
frequently  used  as  a  measure  of  the  charge  separation  in  the  double  layer. 


III.  DENSITY  FUNCTIONAL  FREE  ENERGY  THEORY 
General  formalism 

We  start  our  study  of  the  double  layer  problem  with  a  discussion  of  the  grand  canonical 
density  functional  formalism  for  a  mixture  of  ionic  species  in  an  external  field.  In  this  work 
we  adopt  the  general  approach  due  to  Morita  and  Hiroike  [40],  De  Dominicis  [41],  Stillinger 
and  Buff  [42]  and  to  Lebowitz  and  Percus  [43],  and  used  later  by  many  investigators  [20- 
22,44].  The  particle-particle  direct  correlation  functions  of  a  non-uniform  fluid  which 
appear  in  the  formalism  allow  us  to  write  an  exact  expression  for  the  density  distribution 
functions.  The  same  formalism  has  been  used  by  Fortsmann  and  collaborators  [13]  as  the 
starting  point  of  their  HNC/LMSA  theory  of  the  double  layer.  This  approach,  which  is 
due  to  Mermin  [45]  and  was  employed  by  Hohenberg  and  Kohn  [46]  for  the  inhomogeneous 
electron  gas,  is  naturally  expressed  in  the  language  of  the  grand  canonical  ensamble. 

The  properties  of  an  interface  and  its  coexisting  bulk  fluid  are  determined  by  the 
constancy  of  the  chemical  potentials,  fia,  and  temperature, T,  throughout  the  system.  The 
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free  energy  appropriate  to  the  grand  canonical  ensamble  is  the  grand  potential.  For  a 
mixture  of  particles  of  different  kinds  (a  =  1,...,  c),  the  grand  potential  functional  is 


ft  =  - kTinE ,  (3.1) 

where  E  is  the  grand  pairtition  function  and  k  is  Boltzmann’s  constant.  The  equilibrium 
density  distribution  is  an  unconstrained  minimum  of  the  grand  potential  functional,  ft, 
where 

fi({n})  =  F({n»  -  W  d3rfiana(r).  (3.2) 

a  J 

Here,  F({n})  is  the  Helmholtz  free  energy  functional  of  the  system  and  {n}  denotes  the 
functional  dependence  of  ft  and  F  on  the  particle  densities,  na(r),or  =  l,...,c. 

For  a  mixture  of  charged  particles  of  species  at  absolute  temperature  T  in  the  field  of 
an  externad  potential  vQ(r),  the  grand  potential  functionad  can  be  written  as 

fi(W)  =  rf3rMrXva(r)  ~  Mo) 

°  ^  f  (3-3) 

+  kT^2  j  <f3rna(r)(^n(A’na(r))  -  1]  -  <KM). 


The  second  term  on  the  rhs  of  Eq.(3.3)  is  the  ideal  gas  contribution  to  the  Helmholtz  free 
energy  and  Aa  is  the  thermal  de  Broglie  wavelength  of  patrticles  a.  The  term  -<p  in  the 
same  equation  corresponds  to  the  interpaurticle  interaction  contribution  to  the  free  energy. 

The  grand  potential  functional  ft({n})  is  minimized,  for  fixed  va(r)  and  uas( r,r'), 
when  nQ(r)  talces  its  equilibrium  value.  In  that  case,  ft  corresponds  to  the  equilibrium 
grand  potential  function.  The  functional  - <j> ,  on  the  other  hamd,  cam  be  used  ais  a  gener¬ 
ating  functional  for  n-body  correlation  functions,  in  particular,  from  the  first  functional 
derivative  we  obtain: 

=  c°(r;<n))’  <3-4> 

while  the  second  functional  derivative  of  <t>({n})  defines  the  Ornstein-Zernike  direct  corre¬ 
lation  function 

1  ^(M)  frj.|nn  /« 

kT  6na{r)6n0(r)  c«*(r’r  ’  {n»-  (3,5) 

The  equilibrium  condition  cam  then  be  expressed  ais 

=  kTtn{na{ r)/Ca)  +  Mr)  -  kTca( r;  {n})  =  0,  (3.6) 

where  Ca  =  A”3  exp(/fyt<*)  is  the  fugacity  of  component  a  in  the  mixture  amd  15  =  1/kT. 
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By  functional  integration  between  an  initial  state  n1,  and  a  final  state  n,  it  is  possible 
to  obtain 

t({n})  =  4>({n'})  +  kT  Y  /  d3r[na(r)  -  n^(r)]cQ(r;  {n1}) 

a  J 

+  ’'T'£.J  JMr' ["•(')-<(')!  (3.7) 

a/3 

x  [n5(r') -n^(r')]^  d\  J  d\cQ0(r,  r';  A). 


In  order  to  obtain  this  result,  the  linear  density  path, 


na(r;  A)  =  nj,(r)  +  A[na(r)  -  nj,(r)]. 


was  used  for  the  integration.  The  parameter  A  can  take  values  in  the  interval  0  <  A  <  1. 
Equations  (3.3),  (3.6)  and  (3.7)  can  now  be  employed  to  write  the  following  expression  for 
the  grand  potential  functional: 

fi(M)  =  ft({n’})  +  Y  [  d3rn»(r)[va(r)  ~  uo(r)] 

a  ** 

+  A tfW  d>m.(r)<n(n.(r)/ni(r)) 

a  J 

-  kT  S  /  dMn«(r)  “  n«(r)l  (3.8) 

a  J 

~  kT^2  J  J d3rd3r'[na{ r)  -  nj,(r)][n*(r')  -  nj,(r')J 

x  f  dX  f  d\'ca0(r,  r';  A'). 

Jo  Jo 

When  dealing  with  long-ranged  Coulombic  potentials  it  is  convenient  to  define  a  short- 
range  part  of  the  direct  correlation  function  c%jj( r,  r')  by 

c*0(r,  r')  =  -0qaqsuc(\  r  -  r'  |)  +  cf  J(r,  r').  (3.9) 

The  short-range  correlation  function,  cf^(r,r'),  can  be  further  separated  by  substracting 
from  it  the  hard-sphere  contribution, cW5(r,r').  This  is, 


Aca0(r,  r')  =  cafi(r,  r')  -  c”/( r,  r'). 


(3.10) 
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These  definitions  allow  us  to  rewrite  the  grand  free  energy  functional  as 
ft({n})  =  Q({n‘})  +  j  d3rna(r)[va(r)  -  w£(r)] 

a 

+  kT^T  f  d3rnQ( rj£n(na(r)/njt(r)) 

~kTYl  I <f3r(n«(r) ~  na(r)3 

“  ‘  (3-11) 

+  2  S  J  J  j3r£f3r'[n“(r)  "  no(r)][n^(r')  -  nl^r')\qaq0uc(\  r  -  r'  |) 

a0 

-kTYs  J  J  d3rd3r'[na(r) -n'a(r)][ni3(r') -n'pir') } 

x  [  dX  f  d\'AcaS(T,r';X')  +  AFHS({n}). 

Jo  Jo 

The  last  term  on  the  right  hand  side  of  Eq.(3.11)  represents  the  excess  free  energy  change, 
between  the  initial  state  n\  and  the  final  state  n,  produced  by  the  hard  sphere  interaction 
exclusively  (in  the  presence  of  the  other  interactions). 

Using  the  equilibrium  condition,  Eq.(3.6),  we  can  obtain  the  following  formal  expres¬ 
sion  for  the  equilibrium  density  profiles  nQ( r): 


kT£n(na(r)  /  n'a(r))  =  -«(r)  -  (r))  -  qa(ip( r)  -  0*(r)) 

+  kT  ^2  J  d3r’[n0( r')  -  nj,(r')]  d\Aca0(r,r'-,  A) 

_  SAFHS({n}) 

6na(  r) 


(3-12) 


In  the  derivation  of  Eq(3.12)  use  has  been  made  of  the  definition  of  the  mean  electrostatic 
potential,  Eq.  (2.7). 


Generalized  hard-rod  model 

The  formalism  presented  above  must  be  completed  by  specification  of  a  model  func¬ 
tional  for  the  excess  free  energy  of  a  hard-sphere  mixture.  Generalizing  the  exact  result 
for  an  inhomogeneous  system  of  hard-rods,  Percus  [36]  and  Vanderlick  et  al. [29]  have  de¬ 
fined  a  generic  free  energy  functional  for  the  inhomogeneous  hard  sphere  fluid.  In  three 
dimensions,  the  free  energy  function  of  a  hard  sphere  mixture  is  [29] 


f*“"  =  E  /  ^(rJ-FVaa'fr)}), 


(3.13) 
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where  F0({n(r)})  is  the  excess  free  energy  per  particle  of  a  homogeneous  mixture  of  hard 
spheres  evaluated  at  the  position  r  and  ^"(r)  and  nr(r)  are  coarse  grain  densities.  Each 
one  of  these  densities  is  defined  by  a  weighting  function  of  the  relative  position  to  the  hard 
sphere  center,  and,*in  the  most  general  case,  is  also  a  functional  of  the  density  distribution, 


n£(r)  =  J  d3r'ua(r  -  r';  {n})n0(r'), 

(3.14) 

n;(r)  =  /  dVr^r  -  r';  {n})n«(r'). 

(3.15) 

Use  of  definitions  (3.14)  and  (3.15)  for  a  uniform  fluid  mixture  gives  the  following  normal¬ 
ization  conditions: 


J  d3rua(r  -  r';  {n})  =  J  d3rra( r  -  r';  {n})  =  1.  (3.16) 

In  order  to  establish  a  theory  for  strongly  inhomogeneous  fluids  based  on  Eqs.(3.12)-(3.16), 
a  particular  form  for  the  weighting  functions  v  and  r  must  be  specified.  The  assignment 
of  weighting  functions  generates  different  model  density  functionals.  A  disscusion  of  how 
different  appropriate  selections  of  weighting  functions  generate  several  important  model 
functionals  can  be  found  in  Ref.(30).  Of  particular  importance  for  this  work  are  the  forms 
of  u(r)  and  r(r)  proposed  by  Robledo  and  Varea  [37],  and  by  Fischer  and  Heinbuch  [38]  as 
three  dimensional  generalizations  of  the  hard-rod  model.  This  model,  termed  here  Gener¬ 
alized  Hard  Rod  Model  (GHRM),  is  characterized  by  the  following  weighting  functions: 

►.(r  -  t>)  =  «((<*./ 2)-  |  r  -  r'  |)/(4*(<f„/2)2),  (3.17) 

r„(r  -  r')  =  H((da/2)~  |  r  -  V  |)/(4*(<t,/2)V3),  (3.18) 

where  da  is  the  diameter  of  the  particles  in  the  fluid,  6(r)  is  the  Dirac  delta  function,  and 
H(r)  is  the  Heaviside  step  function: 


H( r)  =  l,  r  >  0, 
=  0,  r  <  0. 


(3.19) 


In  this  model,  the  coarse  grain  density  n£(r)  is  the  average  density  of  species  over  the 
surface  of  a  sphere  of  radius  da/2.  The  coarse-grain  density  fjj(r)  is  the  average  density 
of  species  inside  a  sphere  of  radius  da/2. 

The  functional  in  Eq.  (3.13)  allows  us  to  determine  an  expression  for  the  free  energy 
change  AFW5({n})  appearing  in  Eq.(3.11).  It  is 

AFHS({n})  =  £  J  d3r  J*  rfA  Jj  (n£(r;  \)F0  ((nr(p;  A)})) .  (3.20) 
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Using  this  expression  we  can  finally  rewrite  our  equation  for  the  equilibrium  density  profiles, 
Eq.(3.12),  as 


kT£n(na(r) / n^ir))  =  -«(r)  -  (r))  -  ?a(0( p)  -  0'(*O) 

+  J  d3r'[n0(r')  -  n^(r')]  J  dAAcQ/*(r,  r'\  A) 

4,m"’ 


6na(  r) 


F0({iT(r');A}) 


—  3fa({nr(r';A)})W;(r';A) 
3n^(r')  6n0(r) 


(3-21) 


To  study  a  bulk  fluid  in  equilibrium  with  a  planar  electrode  we  now  identify  the 
initial  state  {n1}  with  the  neutral  bulk  electrolyte.  This  corresponds  to  an  homogeneous 
solution  in  which  no  external  forces  are  present.  Since  we  are  considering  an  infinite  plane 
with  uniform  charge  density  ,  local  densities  vary  only  in  the  direction  x  normal  to  the 
wall.  Additionally,  Eq.(3.21)  requires  the  knowledge  of  inhomogeneous  direct  correlation 
function  in  excess  over  the  hard  sphere, Aca^(r,  r'),  for  all  the  possible  positions  r  and  r' 
across  the  interface.  Since  these  correlation  functions  are  not  known,  we  approximate  the 
function  Aca|9(r,r')  with  the  function  Acajs(l  r  -  r'  |)  of  the  homogeneous  neutral  bulk 
electrolyte  in  equilibrium  with  the  interface,  this  is 

AcQ/3(r,  r';  A)  S  Aca/J(r,  r';  A  =  0)  =  Aco5(|  r  -  r'  |).  (3.22) 

It  is  convenient  to  emphasize  here  that  Aca^(|  r  —  r'  |)  is  a  pair  correlation  function  for 
a  neutral  bulk  electrolyte  whereas  the  interface  is  locally  non-neutral.  With  this  approx¬ 
imation,  and  by  using  Eq.s  (3.14)  and  (3.15)  for  the  GHRM,  Eq.(3.21)  can  be  rewritten, 
for  a  planar  symmetry,  as 


na{x)/na  =  exp{-0qatp(x) 

-  r  ix’vU*  -  i')0F„«a'(i')))  +  0F.( n) 

-£/o  dzn>(z)  «!;(»<—  T*»<*  "  *  >  +  A. TC 

+  YlJ i3r'(n«(x')  -nj]Aco»(|  r-r'  |;{n})}, 


and 


(3.23) 


for  x  <  da/ 2, 


na(x)  =  0,  for  x  >  da/2. 
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For  the  planar  symmetry  the  coarse  grain  densities  n£(x)  and  n£(x)  can 

be  written  as 

*  h£(x)  =  j  vta{x  -  x')nQ(x')dx' , 

(3.24) 

and 

txa(x)  =  J  TXa(x  x  )nQ(x  )dx  , 

(3.25) 

where  the  reduced  weighting  functions  vtQ  and  tzq,  are  defined  by 

t'za(x)  =  J  J  va(r)dydz, 

(3.26) 

Txa{x)  =  J  j  ra(r)dydz. 

(3.27) 

After  integration  over  coordinates  y  and  z ,  we  obtain,  for  the  GHRM, 

^zo(x)  =  H((da/2)-  |  x  | )/da, 

(3.28) 

rza(x)  =  QH((da/2)~  I  X  |)((da/2)2  -  x2)/d3. 

(3.29) 

We  approximate  now  the  excess  free  energy  per  particle  F0  of  the  homogeneous  hard- 
sphere  fluid  by  the  Carnahan-Starling  equation  of  state  [47].  This  equation  was  generalized 
to  a  mixture  of  hard  spheres  of  different  sizes  by  Mansoori  et  al.  [48].  The  expression  for 
the  free  energy  is 


=  -^(1  “  yi  +  V2  +  2/3)  +  (3y2  +  2y3)(l  -  0  1 
+  |(i  -  yi  -  y2  -  ^yaXi  -  0-2  +  (ya  -  i)Mi  -  0 


(3.30) 


where 

£  —  53^01  —  uanai  va  = 

a 

The  variables  y3,  y2,  and  y3  are  defined  as  follows: 

Vi  =  5  T,  +  dff)/(dadff)ll2f 

a/3 

at$  y 

r  i3 


ya  = 


1 


(3.31) 


(3.32) 

(3.33) 

(3.34) 


11 


(3.35) 


and 

Aa/j  =  (vaV0)inan0(da  -  dp)2 /(Zndcdp). 

In  the  last  two  equations,  n  is  the  total  number  density  given  by  J2a  n°- 

In  their  study  of  fluids  confined  between  planar  wails,  Vanderlick  et  al.  [30]  compared 
three  different  approximate  density  functional  free  energy  theories  of  inhomogeneous  fluids 
for  hard  spheres  and  Lennard-Jones  potentials.  Their  study  included  the  GHRM  and  the 
SDA  due  to  Tarazona  [24].  The  results  of  that  study  show  that  whereas  the  GHRM  is 
quantitatively  inferior  to  the  SDA,  it  is  qualitatively  correct.  Since  the  GHRM  captures  the 
qualitative  behavior  of  confined  hard-sphere  and  Lennard-Jones  fluids  and  retains  enough 
mathematical  simplicity  we  are  encouraged  to  apply  it  to  more  complicated  systems. 

MSA  Approximation  for  Aca^(j  r  -  r'  j) 

It  follows  from  Eq.(3.23)  that  our  description  of  the  electrical  planar  interface  is  still 
not  complete  without  a  prescription  for  bulk  phase  direct  correlation  functions  Aca^(| 
r  —  r'  |).  Several  choices  can  be  immediately  invoked  from  bulk  electrolyte  theory.  A 
simple  choice  is  to  use  the  direct  correlation  functions  of  the  MSA.  The  MSA  is  a  relatively 
accurate  approximation  which  generates  analytic  expressions  for  the  direct  correlation 
function  of  several  important  model  potentials  [49].  Waismann  and  Lebowitz  [50]  showed 
that  the  integral  equation  resulting  from  the  Omstein-Zemike  equation,  has  an  analytical 
solution  when  the  MSA  closure  for  the  restricted  primitive  model(RPM)  is  employed.  The 
RPM  is  a  still  simpler  model  in  which  all  ions  have  the  same  size;  da  =  d.  For  the  RPM, 
the  MSA  provide  the  following  expression  for  the  function  Ac0/j(|  r  —  r'  |): 

[(2 B/d)  -  (B/d)*3  -  1/j]  , .  <  d 
=  0  s  >  d, 

where  s  =|  r  —  r'  |  and 

S«[V>  +  1-(1  +  2S?)*]/V)  (3.37) 

and  <t>  =  nod.  The  quantity  is  the  inverse  Debye  screening  length  given  by 

C 

K2D  =  (4n/3/e)Y,naql  (3.38) 

or 

The  solution  given  by  Eq.(3.36)  holds  for  an  arbitrary  number,  c,  of  ionic  species  provided 
global  charge  neutrality  is  maintained, 

C 

TlaQa  =  0.  (3.39) 

a 

The  MSA  direct  correlation  function  is  an  important  piece  in  the  formulation  of  the 
HNC/MSA  theory  of  the  double  layer  [9].  In  that  theory,  the  approximation 

Ca0(r,  r';  A)  3  ca${ r,  r';  A  =  0)  =  caj*(|  r  -  r'  |)  (3.40) 
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is  made;  thus,  the  effect  of  the  external  potential  and  the  inhomogeneities  of  the  interface 
on  c(r,  r')  are  entirely  neglected.  We  believe  that  this  approximation  is  more  severe  than 
the  similar  approximation  of  Eq.(3.22)  used  in  this  work. 


IV.  RESULTS 

In  this  section  we  present  our  results  for  the  density  profiles,  mean  electrostatic  poten¬ 
tial  profile,  and  diffuse  layer  potential  drop  for  solutions  containing  symmetrical  1:1  and 
2:2  electrolytes.  The  results  are  compared  with  existing  MC  data  and,  when  possible,  with 
results  obtained  from  severed  other  approximations.  The  calculations  were  performed  by 
means  of  the  method  of  subdomains,  finite  element  basis  functions,  collocation  weighted 
residuals,  [51]  and  Newton  iteration  with  initialization  chosen  by  parametric  continuation 
[52].  We  choose  quadratic  Lagrange  interpolating  polynomials  as  basis  functions.  This  nu¬ 
merical  technique  was  applied  before  to  the  solution  of  PY,  HNC  [52]  and  MSA  [53]  integral 
equations  for  bulk  simple  fluids  and  was  extended  by  Mier  y  Teran  et  al.  [9]  for  solving 
the  HNC/MSA  integral  equation  for  the  double  layer  problem.  A  detailed  comparative 
discussion  about  the  application  of  this  method  to  the  solution  of  HNC/MSA  equation  for 
the  double  layer  RPM  and  its  efficiency  and  accuracy  can  be  found  in  Ref.(54). 

With  the  algorithm  mentioned  above,  we  reduced  the  set  of  Eqs.(3.23)  to  a  system  of 
algebraic  equations  for  the  values  of  reduced  density  profiles  at  the  nodes:  gai  =  na(xi) /na. 
This  nonlinear  set  of  equations  is  solved  by  Newton’s  method.  The  iterative  process  is 
continued  until  the  Euclidean  norm  of  the  updates  after  iteration  it  +  1  becomes  less  than 
IQ’10: 


.  at  «=1 

where  N  is  the  number  of  nodes  in  the  domain  d/2  <  x  <  R,  and  R  is  the  cutoff  value  for 
he  integrals  in  Eqs.(3-23).  Both  the  number  of  nodes  N  and  the  value  of  R  depend  on 
concentration.  We  used  a  uniform  mesh  in  the  domain  d/2  <  x  <  R. 

Either  the  charge  density,  <r,  or  the  electrostatic  potential  at  the  electrode,  0e,  can 
be  specified  and  the  equations  solved.  At  very  low  charge  densities,  <7,  or  potentials, 
V>«,  we  found  it  convenient  to  use  the  MGC  density  profiles  as  an  initial  guess.  Once  a 
solution  for  certain  values  of  the  parameters  is  found,  initial  estimates  at  other  values 
for  the  parameters  can  be  found  easily  by  a  first-order  continuation  technique.  Typically, 
three  to  five  Newton  iterations  are  needed  to  reach  convergence.  After  convergence  was 
attained,  the  value  of  a  or  the  value  of  xj>e,  depending  on  which  quantity  was  used  as 
the  parameter,  were  computed  using  Eq.(2.8)  or  Eq.(2.10)  respectively.  The  agreement 
with  the  value  of  a  or  t/>e  originally  used  to  solve  the  equations  gives  an  indication  of  the 
accuracy  of  the  numerical  method.  Except  for  the  very  low  concentration  regime,  Eq.(2.8) 


<  10 


-10 


(4.1) 
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or  (2.10)  was  satisfied  to  at  least  five  significant  figures.  In  the  most  dramatic  case  treated; 
1:1  electrolyte  at  0.01  M,  the  dimensionless  Debye  distance,  k^1  ,  becomes  very  small  and 
Eq.(2.8)  or  (2.10)  was  satisfied  to  four  significant  figures  only. 

In  our  calculations  we  have  used  dimensionless  parameters.  We  express  all  lengths 
in  units  of  the  diameter  d.  The  dimensionless  surface  charge  is  a *  =  udP /e,  where  e  is 
the  magnitude  of  the  electronic  charge.  Similarly  the  dimensionless  potential  profile  is 
tp*(x)  =  /?et/>(x).  In  order  to  compare  with  the  MC  data  of  Valleau  and  co-workers,  [6]  we 
fixed  the  value  of  the  plasma  parameter  to  T  =  /3e2/ed  =  1.6809.  This  value  corresponds 
to 

e  =  78.5,  T  =  29 8K  and  d  =  4.25A. 


1:1  electrolytes 

We  solved  Eq(3.23)  for  electrolytes  ranging  from  0.01-2  M  and  surface  charge  ranging 
from  0.05-0.9.  In  Table  I  we  list  the  dimensionless  diffuse  layer  potential,  0eip(O),  where 
t/>(0)  is  the  potential  drop  between  the  point  of  closest  approach  to  the  surface  and  infinity. 
Note  that  the  position  of  the  wall  has  been  shifted  to  x  =  -  d/2.  In  Table  I  we  also  display 
the  MC  results  [6],  those  of  MGC,  BGY  [11],  MPB5  [7]  theories  and  the  OZ/LMBW 
results  obtained  recently  by  Plischke  and  Henderson  [15]  using  the  HNC  closure.  The 
general  agreement  of  our  results  with  the  MC  data  is  quite  good.  A  clearer  comparison  of 
our  results  for  1:1  electrolytes  with  MC  data  is  given  by  Fig.  1  where  we  plot  the  diffuse 
layer  potential  t/>*(0)  as  a  function  of  the  reduced  charge  density  a* .  As  reported  before 
by  other  authors,  [11,15]  in  the  low  concentration  regime,  density  profiles  become  very 
long  ranged  and  special  numerical  difficulties  appear.  We  believe  that  the  discrepancies 
between  our  results  and  the  MC  data  at  c  =  0.01  M  can  be  attributed  ,  at  least  in  part, 
to  this  cause.  The  crosses  shown  in  Fig.  1  are  the  results  of  Plischke  and  Henderson  [15]. 

The  classical  MGC  theory,  which  neglects  the  finite  size  of  the  ions,  predicts  an  interfa¬ 
cial  thickness  which  is  greater  than  that  obtained  by  MC  simulations  for  low  surface  charge 
densities,  and  smaller  than  that  of  the  MC  data  for  large  <7*.  This  phenomenon  is  evident 
at  1  M  concentration.  In  Fig.  2  we  plot  the  diffuse  layer  potential  as  a  function  of  surface 
charge  density  for  c  =  1  M.  In  that  figure  we  also  display  the  MC  data,  and  the  results 
obtained  from  the  approximations  listed  in  Table  I.  The  agreement  of  our  results  with  MC 
is  very  good  and  in  some  cases  of  comparable  accuracy  with  those  obtained  by  Plischke 
and  Henderson  [15]  with  the  OZ/LMBW  and  the  HNC  closure.  The  data  available  for  the 
MPB5  theory  show  an  excellent  agreement  with  MC  data.  Unfortunately  the  data  are  for 
low  values  of  a*  only,  and  because  of  the  secondary  role  played  by  the  excluded  volume 
effects  in  the  MPB5  theory,  it  is  not  expected  that  the  theory  can  be  applied  at  higher 
surface  charges  where  the  size  effects  are  very  important.  The  BGY  theory  of  Caccamo  et 
al.,  Ref.(ll),  wich  is  very  good  for  low  a*  fails  to  predict  the  change  in  curvature  showed 
by  the  MC  results  at  intermediate  charge  densities. 

The  'assical  theory  of  Gouy-Chapman  always  predict  monotonic  variation  for  the 
density  profiles  of  both  coions  and  counterions.  In  contrast,  for  1:1  electrolytes,  the  MC 
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results  of  Valleau  and  collaborators  [6]  for  the  structure  of  the  RPM  double  layer  exhibit 
interesting  layering  effects  for  high  surface  charges.  In  Fig.  3  we  present  a  comparison  of 
our  results  for  the  density  profiles  of  a  double  layer  for  a  bulk  density  of  c  =  0.1  M  and 
<7*  =  0.30,  with  those  corresponding  to  MC  simulation  and  the  MGC  theory.  Our  results 
agree  quite  well  with  the  MC  results.  All  the  profiles  showed  are  monotonic  in  this  case. 
In  Fig.  4  we  plot  the  mean  electrostatic  profile  which  correspond  to  the  same  condition 
presented  in  Fig.  3.  Again  we  obtain  very  good  agreement  with  the  MC  results.  The  MGC 
theory  is  relatively  succesful  in  describing  both  density  profiles  and  mean  electrostatic 
potential  at  c  =  0.1  M  and  a*  =  0.30. 


In  Fig.  5  we  present,  with  solid  lines,  the  counterion  and  coion  density  profiles  obtained 
in  this  work  for  c  =  1 M  and  <7*  =  0.42.  For  comparison,  in  the  same  figure  we  show  the 
MC  and  MGC  results.  Also  shown  in  the  figure  are  the  results  obtained  by  Plischke  and 
Henderson  for  the  OZ/LMBW  theory  with  the  MSA  closure.  It  is  important  to  mention 
at  this  point  that,  for  this  concentration  and  surface  charge,  the  results  of  the  OZ/LMBW 
theory  with  the  MSA  closure  are  in  very  good  agreement  with  those  of  the  same  theory 
when  the  HNC  closure  is  employed  [15].  The  MC  results  clearly  show  the  onset  of  the 
formation  of  a  second  layer  of  counterions  near  x/d  =  1.  Since  the  MGC  theory  is  a 
point  chaxge  theory,  it  does  not  predict  the  layering  phenomenon.  On  the  other  hand,  the 
OZ/LMBW  theory  accurately  follows  the  behavior  of  the  MC  data.  It  is  interesting  to 
see  that  the  density  functional  theory  presented  in  this  paper  is  also  able  to  predict  the 
formation  of  the  second  layer  of  counterions.  However,  the  position  of  the  second  layer  is 
clearly  shifted  towards  the  electrode.  The  theory  also  exaggerates  the  size  of  the  second 
peak.  The  coion  density  profile  predicted  by  the  density  functional  theory  agrees  quite 
well  with  MC  data  and  is  almost  indistinguishable  from  that  of  the  OZ/LMBW  theory. 

At  a  bulk  concentration  of  c  =  1 M  and  a  charge  density  cr*  =  0.7,  a  second  layer 
of  counterions  is  clearly  formed.  In  Fig.  6  we  compare  the  density  profiles  predicted  by 
MGC,  OZ/LMBW,  and  density  functional  theories  with  the  MC  results  for  that  conditions. 
Again  the  MGC  theory  predicts  monotonic  profiles  while  the  OZ/LMBW  theory  very  well 
predicts  both  the  position  and  magnitud  of  the  second  layer.  Again  the  density  functional 
theory  overemphasize  the  value  of  the  density  of  the  second  layer.  As  in  Fig.  5,  the 
position  of  this  layer  is  shifted  towards  the  electrode.  This  can  be  a  consequence  of  the 
way  in  which  the  GHRM  takes  into  account  the  hard-core  effects.  The  GHRM  predicts,  for 
a  hard-sphere  fluid  near  a  hard  wall,  a  density  profile  with  a  second  peak  shifted  towards 
the  wall  when  compared  with  simulation  results.  See  Ref.  (30). 


The  mean  electrostatic  potential  profile  wich  corresponds  to  the  last  conditions  pre¬ 
sented  is  shown  in  Fig.  7.  The  agreement  between  the  density  functional  theory  and  the 
MC  results  is  very  good.  The  density  functioned  theory  is  able  to  predict  the  presence 
of  a  very  shallow  minimun  in  this  function.  A  similar  minimum  is  present  in  the  MC 
results.  The  mean  electrostatic  potential  function  is  not  very  sensitive  to  the  details  in  the 
structure  of  the  double  layer. 
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2:2  electrolytes 

We  have  computed  results  of  our  density  functional  theory  for  two  concentrations: 
0.05  and  0.5  M.  In  the  lower  part  of  Table  I  we  display  some  results  of  this  work  for  the 
diffuse  layer  potential  and  compare  with  those  of  the  theories  previously  mentioned.  We 
find  reasonable  agreement  with  MC  data. 

In  the  case  of  divalent  electrolytes,  the  MC  results  show  the  interesting  phenomenon 
of  charge  inversion.  This  phenomena,  which  is  a  result  of  both  hard-core  and  electrostatic 
interactions,  consists  in  the  formation  of  a  second  layer  of  coions  next  to  the  first  layer 
of  counterions.  In  Fig.  8  we  plot  density  profiles  for  a  double  layer  at  c  =  0.5A/  and 
cr*  =  0.1704.  Nearly  all  the  counterion  charge  is  concentrated  into  a  thin  layer  next  to  the 
wall.  The  response  of  the  system  to  this  dipole  layer  is  the  formation  of  a  layer  of  coions 
within  x  =  d  and  x  =  2d  approximately.  As  can  be  seen  in  Fig.  8,  the  density  functional 
theory  is  predicting  the  charge  inversion  phenomenon.  We  obtained  a  counterion  density 
profile  in  very  reasonable  agreement  with  MC  data.  On  the  other  hand,  our  theory  tends 
to  underestimate  the  magnitude  of  the  maximun  in  the  coion  profile.  As  expected,  the 
MGC  theory  totally  ignores  the  charge  inversionn 

In  Fig.  9  we  show  the  mean  electrostatic  potential  profile  for  the  same  conditions 
presented  in  the  previous  figure.  The  MC  simulations  result  in  a  potential  profile  which 
is  oscillatory  with  a  minimum  of  about  -0.2  just  beyond  one  diameter  from  the  wall.  Our 
density  functional  theory  predicts  the  oscillatory  behavior  and  is  in  very  good  agreement 
with  the  MC  data. 


V.  SUMMARY 

We  presented  a  nonlocal  free  energy  density  functional  theory  for  the  electrical  double 
layer.  Within  the  frame  of  the  grand  canonical  formalism,  we  construct  a  free  energy 
functional  of  the  density  distribution.  We  then  separate  the  short  ranged  part  of  the 
inhomogeneous  direct  correlation  functions  which  appear  in  the  formalism  into  a  hard- 
sphere  term  and  a  residual  term.  The  residual  term  contains  the  correlations  arising  from 
the  Coulombic  interactions  between  particles  in  the  fluid.  The  hard-sphere  part  of  the 
free  energy  functional  is  then  approximated  by  a  generic  functional  proposed  by  Percus 
[36]  as  a  three  dimensional  generalization  of  an  inhomogeneous  hard-rod  system.  We  used 
its  extension  to  mixtures  due  to  Vanderlick  et  a/.  [29].  This  generic  nonlocal  functional 
requires  the  specification  of  two  coarse-  grain  densities.  In  this  work  we  choose  to  use 
the  weighting  functions  proposed  by  Robledo  and  Varea  [37]  and  by  Fischer  and  Heinbuch 
[38]  to  generate  a  GHRM  functional  for  the  free  energy  of  an  inhomogeneous  hard-sphere 
system.  In  our  calculations  we  approximate  the  free  energy  of  a  bulk  hard-  sphere  mixture 
with  the  Carnahan-Starling  [47]  expression.  The  residual  inhomogeneous  direct  correlation 
functions  are  approximated  with  those  corresponding  to  the  neutral  bulk  electrolyte  which 
is  in  equilibrium  with  the  interface.  Use  is  made  of  the  analytical  solutions  of  the  MSA 
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[47], 


The  GHRM  free  energy  density  functional  theory  described  in  Sec. Ill  correctly  de¬ 
scribes  the  physical  features  presented  by  the  MC  simulations  for  1:1  and  2:2  RPM  elec¬ 
trolytes.  For  1:1  electrolytes  the  theory  predicts  the  layering  of  counterions  which  occurs 
when  the  charge  of  the  electrode  is  increased.  Although  the  theory  exaggerates  the  mag¬ 
nitude  of  the  counterion  layering,  predicts  a  diffuse  layer  potential  which  is  in  very  good 
agreement  with  the  MC  data.  For  2:2  electrolytes,  the  theory  predicts  the  charge  inversion 
phenomenon  and  values  of  the  diffuse  layer  potential  which  axe  in  good  agreement  with 
MC  results.  In  general,  there  are  small  quantitative  rather  than  qualitative  differences 
between  the  MC  results  for  the  density  profiles  and  mean  electrostatic  potential  and  those 
obtained  in  this  work. 

Even  when  calculations  with  the  GHRM  density  functional  theory  are  relatively  sim¬ 
ple,  the  theory  competes  in  accuracy  with  the  more  sophisticated  OZ/LMBW  theory  [15]. 
Because  of  its  simplicity,  the  GHRM  theory  requires  only  a  small  fraction  of  the  computing 
time  used  to  solve  the  OZ/LMBW  theory.  For  the  1:1  electrolyte  at  c  =  1  Af,  our  code 
requires  only  240  s  of  a  Cray2  CPU  time  to  calculate  solutions  and  mean  electrostatic 
potential  profiles  at  15  different  values  of  reduced  charge  density, <7*,  with  a  uniform  mesh 
of  241  nodes,  for  example. 

From  a  density  functional  formulation  similar  to  that  presented  here,  Fortsmann  and 
collaborators  [13]  used  in  the  interface  ion- ion  correlation  functions  of  homogeneous  elec¬ 
trolytes  with  non-neutral  compositions.  Instead,  in  our  work  we  are  employing  a  nonlocal 
GHRM  functional  for  the  hard  sphere  part  of  the  free  energy  and  neutral  bulk  electrolyte 
correlation  functions  for  the  residual  electrostatic  part.  Use  of  non-neutral  composition 
residual  electrostatic  correlations  is  left  for  future  work. 

The  results  of  the  GHRM  for  a  hard-sphere  system  near  a  hard- wall,  reported  in 
Ref.(30)  show  a  poor  quantitative  agreement  with  MC  results.  Since  a  very  good  agreement 
between  the  functional  theory  and  the  MC  results  for  the  planar  double  layer  is  reported 
in  Sec.  IV,  one  can  naturally  ask  if  a  fortuitous  cancellation  of  error  is  occurring  when  we 
combine  the  hard-sphere  free  energy  functional  with  the  MSA  solutions  for  the  electrostatic 
part  of  direct  correlation  functions.  The  answer  to  this  question  probably  can  be  given  by 
solving  the  theory  for  a  more  accurate  functional  for  the  haxd-sphere  contribution  to  the 
free  energy.  The  SDA  of  Tarazona  [24]  seems  to  be  a  good  option  for  this  purpose.  We 
hope  to  contribute  to  the  solution  of  this  question  in  the  near  future. 
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FIGURE  CAPTIONS 

Fig.  1  Reduced  diffuse  layer  potential,  0eip(O),  as  a  function  of  the  charge  density,  ad? /e, 
for  1:1  electrolytes.  Solid  lines  represent  the  results  of  the  GHRM  functional 
density  theory  presented  here  for  0.01  M,  0.1  M,  1  M  and  2  M.  Open  circles, 
solid  squares,  solid  circles  and  open  squares  are  the  corresponding  MC  results. 
The  crosses  (x)  are  the  results  of  the  OZ/LMBW  theory  with  the  HNC  closure, 
Ref.(15). 

Fig.  2  Reduced  diffuse  layer  potential,  (3eip( 0),  as  a  function  of  charge  density,  ad? It, 
for  1  M,  1:1  electrolytes.  The  solid  line  represents  results  of  the  functional  density 
theory.  Solid  circles  are  the  MC  results  of  Valleau  and  collaborators,  Ref.  (6). 
Solid  squares  correspond  to  the  MPB5  theory  [7],  open  circles  to  the  BGY  theory 
[11],  and  open  squares  to  the  OZ/LMBW  theory  with  the  HNC  closure,  Ref.  (15). 

Fig.  3  Reduced  density  profiles,  n(x)/n,  for  a  1:1  electrolyte  at  c  =  0.1M  and  a *  =  0.30. 
The  dots  are  the  MC  results.  The  dashed  lines  correspond  to  the  MGC  theory 
and  the  solid  lines  to  this  work.  Note  that  the  wall  is  at  x  =  —d/2. 

Fig.  4.  Reduced  mean  electrostatic  potential  profile  for  a  1:1  electrolyte  at  c  =  0.1  M  and 
<7*  =  0.3.  All  symbols  as  in  Fig.  3. 

Fig.  5  Reduced  density  profiles,  n(x)/n,  for  a  1:1  electrolyte  at  c  =  lM,am  =  0.42.  The 
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dots  are  the  MC  results.  The  dashed  lines  correspond  to  the  MGC  theory,  the 
dot  dashed  lines  to  the  OZ/LMBW  theory  with  the  MSA  closure,  and  the  solid 
lines  to  thjs  work. 

Fig.  6  Reduced  density  profiles,  n(x)/n,  for  a  1:1  electrolyte  at  c  =  1 M  and  <7*  =  0.7. 
All  symbols  as  in  Fig.  o. 

Fig.  7  Reduced  mean  electrostatic  potential  profile  for  a  1:1  electrolyte  at  c  =  LV/  and 
<7*  =  0.7.  The  dots  are  the  MC  results.  The  dashed  lines  correspond  to  the  MGC 
theory  and  the  solid  line  to  this  work. 

Fig.  8  Reduced  density  profiels,  n(x)/n,  for  a  2:2  electrolyte  at  c  =  0.5M  and  <7*  = 
0.1704.  The  dots  are  the  MC  results.  The  dashed  lines  correspond  to  the  MGC 
theory  and  the  solid  lines  to  this  work. 

Fig,  9  Reduced  mean  electrostatic  potential  profile  for  a  2:2  electrolyte  at  c  =  0.5Af  and 
<7*.=  0.1704.  The  dots  are  the  MC  results.  The  dashed  line  corresponds  to  the 
MGC  theory  and  the  solid  line  to  this  work. 
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TABLE  I.  Diffuse  layer  potential 


c 

<7* 

MGC 

MCa 

BGY*> 

MPB5C 

PHd 

This  work 

1  :  1  electrolytes 

0.01M 

0.10 

5.44 

5.05(0.05) 

— 

5.08 

4.56 

5.26 

0.1M 

0.30 

5.34 

4.63(0.03) 

5.0 

4.74 

4.37 

4.76 

1M 

0.10 

1.4 

1.09(0.06) 

1.055 

1.03 

1.06 

1.03 

0.25 

2.79 

2.13(0.05) 

2.31 

2.10 

2.22 

2.18 

0.42 

3.74 

3.08(0.1) 

3.46 

3.02 

3.23 

3.23 

0.55 

4.26 

4.15(0.15) 

4.21 

- 

4.22 

4.12 

* 

0.60 

4.43 

4.38(0.11) 

4.48 

- 

4.68 

4.52 

0.70 

4.74 

5.71(0.14) 

5.02 

- 

5.76 

5.41 

2M 

0.396 

2.99 

2.29(0.09) 

2.303 

- 

2.29 

2.19 

2  :  2  electrolytes 

0.05M 

0.20 

2.61 

1.33(0.02) 

1.81 

1.36 

1.18 

1.59 

0.5M 

0.1704 

1.36 

0.63(0.04) 

0.64 

0.537 

0.69 

0.57 

G.  M.  Torrie  and  J  P.  Valleau,  Ref.  (6);  1980  and  1982.  Statistical  uncertainty  is  shown  in  parenthesis. 
6  C.  Caccamo,  G.  Pizzimenti,  and  L.  Blum,  Ref.  (11);  1986. 
c  C.  W.  Outhwaite  and  L.B.  Bhuiyan,  Ref.  (7);  1986. 

**  M.  Plischke  and  D.  Henderson,  Ref.  (15);  1988. 
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